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Abstract 

We follow up on our previous works [111 112] which presented a possible approach for deriving 
symplectic schemes for a certain class of highly oscillatory Hamiltonian systems. The approach 
considers the Hamilton- Jacobi form of the equations of motion, formally homogenizes it and 
infers an appropriate symplectic integrator for the original system. In [111 I12| the case of a 
system exhibiting a single constant fast frequency was considered. The present work successfully 
extends the approach to systems that have either one varying fast frequency or several constant 
frequencies. Some related issues are also examined. 
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1 Introduction 



The purpose of this work is to describe an approach for constructing symplectic numerical integrators 
for systems within the following class of highly oscillatory Hamiltonian systems: 

dq dH e dp dH e 

dI = ^p- {q > P) > -dt=-^ q - {q ^ (1) 

with 



H e {qi,q2,Pi,P2) = —^- + ^- + V{qi,q 2 ) + , (2) 



Pi Pi , pIpi , v( s g 2 r »(gi ) 2 92 

where £ <C 1 is a small parameter, where q — (91,92) G K s x and p = (pi,P2) G M. 3 x ^ (the 
upper indices s and / designate the slow and fast variables respectively), and where the interaction 
potential V and the fast frequency factor il(gi) are all independent of e. We assume that f2(gi) is a 
symmetric matrix whose eigenvalues are positive and bounded away from zero and that V is bounded 
from below. We also assume that the initial conditions for ([1]) scale with e so that the initial energy is 
bounded independently of e. Since the small parameter e imparts fast frequencies on the model, direct 
numerical simulation of ([T]) to times of O(l) and beyond is computationally expensive. This type of 
Hamiltonian system has already been studied in many works (2j El El EH HI] • See [4] [12] for a short 
review of the literature, and [TOJ Chap. XIII and XIV] for a general overview. This paper continues 
the work initiated in [111 112) by extending the generality of examples treated in the Hamilton- Jacobi 
framework. The major assumption remaining on the form of ([2]) is that the leading order behavior in 
the fast variables is that of a harmonic oscillator 

plP2 gj^(gi) 2 g2 
2 2e 2 

In Section 2] we briefly discuss the use of coordinate transforms for well-behaved, though fully non- 
linear, potentials [TO] Sec. XIV. 3]. The remainder of that section describes the use of our approach 
on a particular Hamiltonian whose fast potential does not take the form of a harmonic oscillator in 
the original coordinates. 

As in our previous works [Hi 112) , our approach to the problem is based on the Hamilton- Jacobi 
formalism. Let S £ (t,q,P) be the solution to 

d t S e =H e (q + d P S e ,P), S e (0,q,P) = 0. (3) 

For all (q,p, t), it is known that the functions (Q(t),P(t)), implicitly defined by 

P = P(t) + -^(t,q,P(t)), Q(t)=q+-£±(t t q,P(t)), (4) 

are solutions to ([1} with initial conditions (q,p). For any approximation S £ of the generating function 
and stepsize h, there is a corresponding symplectic map '■ (q,p) ^ {Q(h),P(h)) defined by the 
implicit relations 

p = P(h) + ^(h,q,P(h)), Q( h )=q+^(h,q,P(h)). (5) 

In the following we construct a function S s that approximates the solution of ([3]) for small t and e. 
As long as we solve ([5]) exactly, this results in a symplectic numerical scheme (see [6]). Being able to 
generate a class of symplectic schemes motivates the strategy, which we use in the sequel, of making 
all approximations on the level of the generating function S e and from there solving ([5]) to build the 
numerical scheme. Of course, we cannot expect to solve these equations exactly, but we do so to high 
precision and observe good behavior with respect to preserving invariants. 

We work within the parameter regime e <C h <C e 1 / 3 . This yields a computational speed-up in 
comparison to standard algorithms such as velocity Verlet, where the time step must be taken smaller 



2 



than the characteristic time scale of the fast motion, 0(e). However, we note that the energy preser- 
vation property of symplectic schemes is typically proven in the limit h —r for a given Hamiltonian, 
and we are thus working outside the theoretical regime of symplectic schemes. We nonetheless choose 
symplecticity as a goal for our scheme and provide numerical tests of the preservation of energy and 
other invariants. 

In contrast to I12j , where only the case of a system exhibiting one single constant fast frequency 
was considered, the present work successfully extends the approach to systems that have either one 
varying fast frequency, or several constant frequencies. Many of the results presented here have been 
announced in [5]. 

In addition to extending the approach to these two more general settings, we address one issue 
that we left pending in our previous works 112] . This issue is related to the approximation of 
the high order derivatives present in the numerical scheme. One drawback of symplectic integra- 
tors constructed using the Hamilton- Jacobi approach is that they require high order derivatives of 
the Hamiltonian. The question then arises of whether or not these high order derivatives can be 
approximated by the corresponding finite differences without sacrificing the other advantages of the 
approach (symplecticity, accuracy, . . . ) . We demonstrate here that it is indeed possible to use such 
finite difference approximations and keep all the features of the algorithm up to a chosen order. 

The present article is organized as follows. In Section^ we consider a Hamiltonian of the form 
it ( \ Pill , PzP 2 li;/ \ i ni ^2 12 12 

H s {qi,q2,pi,P2) = —^- + —^- + V{qi,q 2 )+0,{qi) (6) 

where fl(qi) is a scalar fast frequency that depends on the slow variables. As was already mentioned 
for the constant frequency case in |12) . the best option to implement our approach is to precondition 
the fast motion using a change of variables. This is to be performed prior to writing the Hamilton- 
Jacobi equation associated to the Hamiltonian dynamics. However, the dependency of fi upon the 
slow variables introduces substantial new difficulties in this preconditioning step as compared to |12) , 
which will be described in details below. 

The algorithm that we obtain following this strategy has been introduced in [5] . In this paper we 
provide a comprehensive discussion of its derivation and properties as well as numerical tests demon- 
strating that the algorithm has favorable error performance in terms of resonances. The algorithm's 
computational efficiency is comparable, although slightly lower due to implicitness, to another estab- 
lished approach for small e (such as 10 -4 ) and does not break down as e — > 0. We further show that 
the exchange of actions (energy divided by frequency) among the fast degrees of freedom is captured 
well by our algorithm. The algorithm presented has many possible variants based on the choice of 
approximation of the generating function in ([3]). Among the possible variants, are those that reduce 
to the algorithms given in |12j whenever fi is constant. 

In Section [3j we next consider the case of a non-scalar constant fast frequency fi, 

tj , \ Pi Pi , P2P2 . s . q 2 ^ 2 q 2 

H e (q 1 ,q2,p 1 ,p 2 ) = ~y~ + ~Y~ + V q2 > + 2e 2 ' ' ' 

We consider first the case when the fast frequencies present in f2 are non- resonant, and next address 
the case when some frequencies are resonant. These terms will be defined in Section [3J We will show 
the efficiency of the algorithm obtained on a classical test case. 

We conclude our work by studying an example system composed of a single point attached by a 
spring to a pivot in the plane (see Section |4]) . Such a system cannot initially be written as a system 
with Hamiltonian of the form After a change of coordinates, we transform the Hamiltonian into 
the form 

H s (qi,q2,Pi,P2) = g + — ~ + v {1ul2) + ^ (8) 

with a non-constant mass matrix M(q 2 ). Form ([5]) is now compatible with @, up to the fact that 
the mass matrix M(q 2 ) depends on the fast position q 2 . We show that our general approach again 
applies, yielding an algorithm with very good numerical performance. 
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Notation Before proceeding, we briefly fix the following notation. For vectors u, v £ R™, we define 
the dot product u ■ v £ R by 



T 

u ■ v — u v 



where we also use as shorthand u 2 = u ■ u. For u € R m and v £ R™, we define the tensor product 

u®v £ R mx ™ by 

(u <g) v)ij — UiVj, for 1 < i < m, 1 < j < n. 
For matrices A,B£ R mx ™, we define A : B e M by 

m n 
i=l = 1 

We employ two different derivative notations. For functions S{q\ 1 P{) and V(qi,q 2 ) with <?i,-Pi £ R s 
and 52 G R^, we denote the entries of q\ by q\ — (<?i,i, . . . , qi, s ), and we write d qi S and ViV^ to denote 
the vectors 

(d<nS)j = jr—, (ViV{q 1 ,q 2 )) j = 7 r^(qi,q2) 1 < j < a. 

2 The case of a scalar non-constant frequency 

In this section, we consider the case © of a Hamiltonian with a scalar non-constant fast frequency 
Q(qi). Slightly changing the notation in comparison to ©, the Hamiltonian reads 

He{qi,q2,Pi,P2) = —^- + —^- + V{qi,q 2 ) + Q{qi) (9) 

We first precondition the equations with a change of variables that takes into account the most 
oscillatory component of the fast solution. Following [TU1 page 555], we introduce the symplectic 
change of variables 



_ Vin(gi) . r . y/n(gi) . y/E . 

qi =qi, pi=pi- or) ,. x q 2 p 2 , qi = — ^— <?2, P2 = — ^==P2, (10) 



which transforms @ into 



Hi{qi,q 2 ,Pi,P2) = ~ PH or)/ \ +^ (P2 P2 + ?2 92j + k gi, g 2 



We neglect part of the slow momentum term and rewrite the slow potential to arrive at 
H 2 (q%,q 2 ,pi,p 2 ) = \p1P1 + ^^-{pIpi +q 2 q 2 ) + V{qx, y/eq 2 ), 



where 

V{qi,q 2 ) = V [qi,—J=q 2 



The system of equations corresponding to H 2 is 



5- t • - Vififgi 1 , T , 
q\ =Pu Pi = -ViK(gi, y/eq 2 ) — — (p 2 p 2 + q 2 q 2 ) 

g 2 = j>2, P2 = -ye v 2 K(gi, veg 2 ) 92 
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In the spirit of [12] . we introduce the fast time 



1 '* 



a(t) = - / Sl(qi{s))ds (11) 
£ Jo 

and consider the time-dependent change of variables (92,^2) >-> (2?, y) defined by 

c/ 2 = xcos(o-(t)) + y sin(tr(t)), p 2 = -xsin(a(t)) + ycos(a(t)), 

which captures the highest frequency component of the fast variables. The transformation is motivated 
by the fact that if V were independent of (72, then x and y would be constant. This gives the dynamics 

<7i = p — > Pi = — 5 o (x x + y y), 

dpi dqx 2e , . 

ac = -^r, 2/ = — — 
ay 

where 

H 3 (qi,x,pi,y,t) = -p\ + V [q x , y/e(x cos cr(t) + ysino-(i))] . 

Note that the dynamics on (51 ,x,pi,y) includes a memory term due to the dependence of er on the 
history of qi, see (|lip . On the other hand, in the case of a constant frequency, fi, there is no memory 
term since a does not depend on e/i, and (I12|) is the Hamiltonian dynamics associated to H 3 . In the 
present case, rather than operate on a system with memory, we consider a as an additional variable 

and add the corresponding conjugate variable a — — (x T x + y y). We then have the dynamics 



9i =Pi, 

Pi = — ViV(gi, i/e(a?cosCT + y sine)) : — <i. 



Vin( g i) 



e 

x = ^/esma V2V(qi 1 \/e(xcosa + ysma)), 
y = — \/ecosa V2V(q\, \fe{x cos a + y sin a)), 

e 

5 = \/e (xsincr — y cos<r) T V2V A (qi, v^xcoscr -I- y sincr)), 
which is a Hamiltonian dynamics associated with the energy 

H4,(qi,x,a,pi,y,a) = \p[p\ + V(c?i, */e(xcos(t + y sin a)) + a ^ g1 ^ . (13) 
2 e 

In the sequel, we take the above Hamiltonian H4 as a starting point for our manipulations. We 
hence write the Hamilton- Jacobi equation associated to H4 and then rescale the variables a, x, and y, 

a x y 

a = — . x = —=, and y = —=, 

e V e V £ 

so that (c/i, x, P, Y, a) are of the same order. We choose not to rescale E despite the fact that it is 
0(t/e) at time t because it plays the role of a fast time in the original dynamics. The Hamilton- Jacobi 
equation becomes 



d t S e (t, qi,x,T,,Pi,Y,a) = #4 

1 



gi + d Pl S e , sfs ( x + -dvS £ ) , S, Pi, y/eY, ea - 



= ^Pi +V(qi+d Px S e ,(ex + d Y S s )cos'S + eYsm'E) ] ^ 



S e (0,gi,x,E,i\,F,o) = 0. 
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The equations corresponding to Q are 

Qi = <7i + d Pl S E , 
P 1 =p 1 - d qi S e , 
where all derivatives of S £ are evaluated at (t, qi, x, E, Pi, Y, a). 
2.1 Expanding in e 

We wish to expand the generating function as a series in powers of e. Recall that there are two 
balancing factors in choosing an ansatz. Since the expanded Hamilton- Jacobi equation of order e k 
will involve the generating function up to order e k+1 , we need to make sufficiently strong assumptions 
to be able to separate out the 0(e k+1 ) terms from the lower order terms into independent equations 
and thereby close the hierarchy. On the other hand, we also need to have a sufficiently general form 
in order to satisfy the resulting equations. We make the two-scale ansatz 

S £ (t,q x ,x, Y,,Pi,Y,a) = S (t, qy , P\ , a) + ^e k S k ^t, E- ir(t, qi, Pi, a), qi, x, £, Pi, Y, a^j (16) 

where we assume S k (t, /?, q±, x, 7, Pi, Y, a) is 2n periodic in 7 and choose the fast time r to satisfy 

d t T=(ViV(qi+dp 1 So,0)+aVin(qi+dp 1 S o ))-dp 1 T + n(qi+d Pl S o ), 
r(Q,q 1 ,P 1 ,a) = 0. 

The periodicity assumption on 7 is needed to separate the orders in the resulting hierarchy of 
equations and is justified both by the fact that equation (|14|) is 27r-periodic in E and that the generating 
function of a harmonic oscillator is periodic in time. The simplifying choice of 5*o independent of E is 
consistent with (|15l) and the fact that A is 0(1) over long times. Similar arguments justify the choice 
that So is a function of (t,qi, Pi,a) only. Note that, as is common in homogenization expansions, 
we have introduced a fast time t/s in Sk, which is considered to be independent from the slow time 
t. Because of the non-constant frequency fl(qi), we choose for r to depend on the slow variables 
(t, qi, Pi, a). What is unusual is that we introduced the fast time as part of a transport term E — t/s. 
We motivate this by briefly considering the Hamiltonian H(a, a) = — ecos(cr) — a/e, a highly simplified 
version of H4 (see (1131) ). It is periodic in a and keeps the conjugate variable a, with frequency = 1. 
In this case, the generating function S e (t, a, E) = — at/e — e 2 [sin E — sin(E — t/e)} satisfies the Hamilton- 
Jacobi equation dtS — H(a + dsS, E) with initial condition 5(0, a, E) = 0. Inspired by the transport 
structure in sin(E — t/e), we choose the independent variables E and E — t/e for Sk, rather than E 
and t/s. 

In view of (|16j). we then have 

dtS e = EZo £k (d t s k -d t T d p s k+ i), dyS e = EZi £kd Y s k, 

d Pl Se = EZo eh ( d Pi S k-d Pl T d p S k+ i), d^S £ = EZi^i^Sk + dySk). 
From the initial condition S E (0, q±, x, E, Pi, Y, a) = and the ansatz, we have 

So (0,<jfi, Pi, a) =0 and S k {0, 7, q u x, 7, Pi, Y, a) = for 1 < k < 00, (18) 
where we note the repetition of arguments in S k results from (|16[) and the fact that r(0, qi, Pi, a) = 0. 



X = x + -dvS e . 

£ 

Y = y- -d x S e> 

e 



E = a + -d a S £ , 

e 

A = a ds^e, 

e 



(15) 



G 



2.1.1 Order e° and e 

Inserting the ansatz (JTH]) into (IT?)) and expanding in terms of e, the O(e ) and 0(e 1 ) equations are 

0(1) : d t S Q - d t r dpSx = \pfPi + V+{a~ d Si - fySifi, 

0(e) : fySi - d t r d p S 2 = V 1 V- (d Pl Si - d Pl Td p S 2 ) 
+ \7 2 V ■ ((x + (9yS'i)cos7 + Fsin7) 

+ (a - d S! - a^v^ • (a Pl 5i - o Pi t d p s 2 ) - (d s 2 + d 7 s 2 )Q, 



where 



V = V(qi+d Pl S -d Pl T dpSi,0) and ft = fi( ?i + 8 Pl S - d Pl T dpS x ), 



As V and SI depend on Si, we must proceed carefully in closing the equations. In fact, we will show 
Si = and close the equations on So by expanding the unknowns So and Si in powers of t. 

Lemma 2.1. The 0(1) and 0(e) equations imply 

dtSo = \p?Pi + V(qi + d Pl S o ,0) + oOfa + d Pl S ), (19) 

Si = 0. (20) 

Proof. We expand the functions So, Si and r in powers of t : 

t 2 t 2 t 2 
So = So,o + tS ,i + ySo,2 H , Si = Si,o + tSi.i + ySi :2 H , r = t + in + -^-t 2 H 

where from the initial conditions of (TIT)) and (fig)) , we have that So,o = tq = 0. The O(£°i ) equation 
is 

S ,i - ndpSi,o = \pfPi + V(qi,0) + (a- d p S h0 - OySifiMqx). 
We wish to separate the So and Si terms. We first observe that (fTT)) gives t x — ft(qi) which leaves 

So,i = ^P?Pi + V(qi,0) + (a -5 7 5i, )fi(gi). 

We integrate both sides of the equation with respect to 7 between and 2tt and note that, by the 

periodicity assumption in the ansatz, +• d^Si t o d~f — 0. Since So is independent of 7, this leaves 

Jo 

So,! = ^P?Pi+V(q u Q) + an(qi), 

which is (fT9")l up to 0(t) error. This implies cLS^o = which, combined with the initial condition on 
S%, implies S^o = 0, so that ([20)) holds up to 0(t) error. 

From here we proceed by induction. Suppose that (fT9)) and (|20l) hold up to 0(t k ). We note 
from (jnj that 

d Pl r = o(; 2 ), 

and thus c^t dpSi — 0(t k+2 ). We then infer from the 0(1) equation that 

d t So - ^idpSi, k = \pIPi + V(qi + d Pl S , 0) + afi(gi + 9 Pl S ) 



fc! 



(a 7 S 1 , fe + 9, 9 S 1 , fe )fi( ?1 ) + 0(^ +1 ). 



Upon canceling out the two dpSi^ terms using ti = f2((?i) and integrating with respect to 7, we find 
that (0 holds up to 0(£ fc+1 ) and that <9 7 Si, fe = 0. 
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Because Si vanishes up to an error of 0(t k ), we infer from pT|) that 

d t r - (ViV + aVifi) ■ d Pl r - = 0(t fe + 2 ). 

We first use this in the O(e) equation to cancel the dpS 2 terms. We then neglect all terms of 0{t k ) 
to find 



t k-i 



(k 



Yy Si, k = V 2 V(qi + d Pl S , 0) • (xcos 7 + Y sin 7) - Sl( qi + d Pl S )d y S 2 + 0{t k ). 



Integrating this equation with respect to 7 gives Si t k = 0. Thus, we have shown that (I20[) holds up to 
0(t k+1 ). Consequently, we have shown by induction that (fT9|) and ([20]) hold to all orders. □ 

In view of (|20p . we have the following simplified expressions for V and Q: 

V = V(q 1 +d Pl S ,0) and Q = Q( qi + d Pl S ). 

In the proof to Lemma |2~T1 we have found 6*0,1 = \Pi P\ + V(qi,0) + afl(qi). Using (fl9|). we likewise 
find 5*0,2 = ViF(gi,0) • Pi + aVif2(gi) • P\. For our generating function, we keep these two lowest 
order terms. We combine them, neglecting only 0(t 3 ) terms, to give 



S Q = t 



-/V/'i i-r( r : -/Vii) +no{ qi + i Pl 



(21) 



We likewise keep the first two orders in the fast time r. From (|17l) . we have To = 0, t\ = £l(qi), and 
t 2 = Wft(qi) ■ Pi. We combine these, neglecting only 0(t 3 ) terms, and take the approximation 

t = m(qi + ^Pi). (22) 



2 

Inserting 5*1 = into the 0(e) equation and using equation (|17p leaves 

d 1 S 2 = ViV- (a; cos 7 + y sin 7). (23) 

2.1.2 Order e 2 

The 0(e 2 ) equation of (fj?]) is 



d t S 2 - d t r 0pS 3 = ~(V U V + aV u O) : (d Pl r ® 9p 1 t)(9 /3 5 2 ) 2 
+ (ViF + aVifi) • (d Pl S 2 - 9 Pi t fyS 3 ) 



Vi 2 U: (^62 «9 Pl T®£) + -V 22 £: W + V 2 ^-drS2COS7 



(24) 



+ (8^2 + d p s 2 ) Vifi • (a Pl r fys 2 ) - (d 7 s 3 + dpS 3 ) n 

where i = a; cos 7 + Y sin 7. Since our approximation r given in (1221) satisfies (|17p up to 0(t 2 ), we can 
cancel the terms involving dp S 3 up to 0(t 2 ), which leaves 

d t S 2 = ^(VuV + aVuh) : {d P J <g> d Pl T){d p S 2 ) 2 + (ViF + aVifi) • d Pl S 2 

- V 12 1> : (9,362 5 Pl f®£) + iv 22 T>:£®£ + V 2 T>-9y52Cos7 
+ (d 7 S 2 + ^5 2 ) Vifi ■ (d Pl r fySa) - d 7 S 3 ^ + 0(t 2 ). 

As before, we integrate with respect to 7 between and 2tt to remove the c^S^ term, leaving only an 
equation in 5*2. This time, however, S 2 is not itself independent of 7. From (1231) we can write 

V 2 V 

S 2 (t,/3,qi,x,~/,Pi,Y,a) = — ^— • (x sin — Y cos 7) + C(i, (3, qi,x, P\,Y, a). 
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The initial condition S^O, 7, q\, x, 7, Pi, Y, a) — implies that 

C{0,/3, qi ,x,I\,Y,a) = - V2 o/ g \ ,0) ' (^sin/3- Fcos/3). 

S%l) 

We integrate ([24]) with respect to 7 between and 2ir, which gives 



(25) 



d t C = V22V 



(x(g> x + Y (g>Y) (V 2 y) 5 



29. 



(Viy + aVi^) -d Pl C 



V 1 Q-d Pl T + |(VuV + aViifi) : {d Pl T®d P ^r) {dpCf + 0{t 2 ) 



where we have used 





p27T p2TT 

j- dyS2 cos 7 d'y = + 
Jo Jo 



JO 

h 



cos 2 7 + dyC cos 7 



d^ 



Si d'y, 
29 ' 



and 



(26) 



2,7V -j^ 1 p27T 

-V22V : I d'y — -V22V : f (x ® x cos 2 7 + a; <E> F cos 7 sin 7 + y <g> Fsin 2 7) d-y 
2 2 J 



1 



V 22 y : (a;®x + y OF). 



We expand 



Cft&gi.s.Pi.y.o) =C , (/9,?i,x,i\,y,o)+tCi()8,gi,s,Pi,y,o) + -C 2 {P,qi,x,P l ,Y,a) + ■ ■ 
and insert this expansion into (|26|) . The O(t) term satisfies 

Cx = Wfe, 0) : (* <8 x + Y ® y) - ^g^ . 
4 2Q(gi) 

We keep the O(t ) term l|25p along with the 0(t) term to approximate S% up to 0(i 2 ) error with 

~ v 2 y(gi + tPi,o) , . v s v 2 y(gi,o) , . „ v R , 

S2 = — — — ; — • {x sin 7 — y cos 7) — — - — • (a; sin p — Y cos p) 



n(gi) 
(v 2 y( gi ,o)) 2 ^ 



-v 22 %o):M, + y 8 F)- 2n(ft) 



(27) 



2.2 Generating function and algorithm 

Combining our previous approximations pip and (|27j> . we approximate the solution to (|14j) by 



1 



2 P^A + y( 51 + ^,0]+ on U + ^i 



+ e 



V 2 y(gi +*Pi,0) 



n(«i + tPi) 

v 2 y(gi,o) 
n(«/i) 



(x sin S — y cos E) 



.r sin ( £ - -0(5! + 2 P i)) ~ rcos ( s _ + 



i (v 2 y) 2 

+^V 22 y(gi, 0) :{x®x + Y ®Y)- t K ^ ' 



(28) 



where we note 



S e (t) = S e + 0(t 6 ) + 0(e 2 t 2 ) + 0{e A ) 



The existence of derivatives of V in the generating function is problematic because it leads to 
higher-order gradient computations in the resulting numerical scheme. In many cases of interest, 
these higher-order gradients are computationally expensive. Furthermore, computing them potentially 
involves extra effort during implementation as most traditional schemes only involve first derivatives of 
V. To avoid this issue in our numerical scheme, we replace all derivatives of V found in the generating 
function with finite differences. We make this replacement before inserting the generating function 
into ([5]) in order to maintain the symplecticity of the scheme. Note that there is no unique way to 
form the finite differences. 

(v 2 y) 2 



In doing this, we choose to exclude 



2ft 



from the 0(e t) term as it is much more computa- 



tionally expensive to implement as a finite difference than the other terms. Rather, we only retain a 
portion of the 0(e 2 t) term that turns out to be sufficient to reproduce the exchange of fast actions 
(see the definition (|33l) and simulation in Fig. [3] below). Note that the only change in the order of 
accuracy is caused by not retaining the portion of the 0{e 2 t) term mentioned above. We choose the 
approximation 



+ 



ipf P a +VU+ |P X , J + ail L + |p a 



V(qi + tPi,e(xsmI] — FcosS)) — V(qi + £Pi,0) 



+ 



t 



Q( qi + tPi) 

e 



V(q u 0) — V [ qi,e [ x sin(S - -fl( qi + 1 - Y cos(E - -f2(<Zi + &)) 

el el 



V(q u ex) - 2V(gi,0) + V(qi,-ex) + V(q u eY) - 2V(q u 0) + V(q u -eY) 



Substituting 5 e .fd into ([5]), we arrive at the following expressions, where we use a = S O (qi + — Pj 
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We first solve for (Pi, Y, E) in the implicit equations 

P±=Pi- h[V 1 V(q 1 + ~i\,0) + aViJlfa + ^Pi)] 

ft. . V2V r (gi,e(x sin a — Y cos cr)) . . . 

- eftVi^fgi + -Pi) v v — r ^ a; cos a + Y sin cr 

2 O(gi) 

ViF (gi + hP x , e(x sin E - Y cos E)) - ViV (gi + hP u 0) 

~ e n( qi + hp 1 ) 

ViO(g x + hPi){V{q x + hP\ , e(x sin £ - y cos E)) - V{q x + hP u 0)) 
+ £ ^(gi + ftP) 2 

ViV(c/i,0) - y 1 V(q 1 ,e(xsina - y cos cr)) 

Vifi(gi)(V(gi , 0) - V( qi , e(x sin a - F cos cr))) 



n( gi )s 



h 



jfriVfaex) - 2ViV(g 1 ,0) + ViV{q u -sx) 
+ ViV(q u eY) - 2V 1 V(q 1 ,0) + ViVfai, -eY) ) , 



V 2 y(gi + hPi,e(xsin£ — FcosE)) . V 2 V(gi,e(x sin cr — Y cos cr)) 

r = y — e — , „ x sin E + e- 



(29) 



fi(?i + /iPi) 



- (v 2 F(gi, ex) - V 2 V(<?i, -ex)) , 



£ = a + ^Q(gi + -P x ). 
We next compute (Qi,X,A) using 



>i = gi + /iPi 
ft 2 



eyViOfe + -Pi) 



V x y(gi + ^Pi,0) + aVififai + ^Pi) 
h V 2 y(gi, e(x sin cr — Fcoscr)) 



0(gi ) 



[x COS cr + Y" sin cr] 



/l£ 



he 



ViV(gi + fePi , e(x sin E - y cos E)) - ViF (gi + ftP x ,0) 
^(gi + ftPi) 

(F(gi + /iPi,e(xsinE-y C osE)) - V( qi + ftPi, 0))Vi^(gi + hPi) 

n( gi + hp 1 ) 2 



V 2 y(gi + /iPi,e(xsinE — y cosE)) VaVYgi, e(x sin cr — Y cos cr)) 
X = x — e — COS E + £ — — ; cos cr 



(30) 



A = 



n(gi + hPi) 

~(y 2 V(qi,eY)-V 2 V(qi,-eY) 

, V 2 F (gi,e(x sin a -Y cos cr)) 
a + £ — — ; |x cos cr + Y sm crj 



O(gi) 

,V 2 V(gi + /iPi,e(xsinE-ycosE)) 
' n(gi + ftPO 



[x cos E + y sin E] , 



where we recall cr = E il ^c/i + —Pi j . The equations (2!)) arc implicit in Z = ( l\ . Y. E) and can 

be written in the form 



Z = z + hF{Z) + sG(Z) + -K(Z) = z + A(Z), 



(31) 
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with z — (pi, y, a) and where K{Z) — (0, 0, fl(qi + t^Pl)) (the dependency of F, G and K on (qi,x, a) 



is not explicitly written). We solve (|31[) for Z with a simple fixed point iteration. In our simulations 
the gradient of A with respect to Z is small since the parameters h and e are small and since in practice 
we work with h 2 /e < 0.4 (see e.g. Fig. 0] below). We thus expect the fixed point algorithm to quickly 
converge. In the test case considered below, this is indeed the case. We terminate the iteration when 
successive iterates differ by a relative factor of 10~ 10 or less, and the algorithm typically converges in 
3 to 8 iterations, depending on the stepsize h and the parameter e. 

The equations CEU]), and O yield Algorithm Q] outlined below. 



Algorithm 1. Initialize: From the initial conditions (gi(0), <?2(0),pi(0),p2(0)), compute the initial 
conditions 

9i (0) = qi (0), Pl (0) = fc(0) - S^^C)' 

2S2(gi(0)) 

Vr>(gi(0)) 1 . 

^(0) = " 92(0), y(0) = - 7 == P2 (0), 

1 ( -77n\- fnl , ^(?l) 2 -T, 



"(0) = 0, O(0) = ^-y [Pi (0)p 2 (0) + (0)&(0) 

Iterate: for rt > 0, 

1. Set (q 1 ,x,a,p 1 ,y,a) = (gf , a;™, <7™,py, y", a n ). 

2. Solve the implicit equations ([2SJ for (Pi, Y, £). 

3. Compute (Qi,X,A) using {3Qj. 

4. Set ( g ?+\z«+V«+\p?+\j/«+ 1 ,a»+ 1 ) = (Q 1; X, E, Pi, Y, A). 

Post-process: From (q^, x , a N ,p^ , y , a ) , return to the original variables: compute first (q 2 1P2) 
using 

= [ x w cos(/) + j/" sm(/)] , p% = yfl [-x N sin( ( r Ar ) + y N cos(a N )] 
and compute next {qf , ,Pi ip2 ) using 



TV AT -JV _ JV , V l^(gj V ) iV JV -JV _ „JV ^JV _ W) W 



As a variant, we will also consider a "no-loop" version of Algorithm [T] Instead of solving (|31j) for 
i?, we perform the following update: 



z* = z + hF(z) + -K{z), 

e 

Z* = z + hF(z*) + sG(z*) + -K(z*), 



h ^ 



e 



and approximate the solution Z = (Pi, Y, S) to (05) (namely, p9]) ) by Z* = (Pf, Y*, £*). Using Z*, 
we next determine (Qi, X, A) in an explicit fashion, using (|30j) . This yields an explicit scheme (we 
first determine z*, next Z* and finally (Qi,X, A)), with a lower computational cost than (|2T)1) - (|3T)1) . 
However, this scheme is not symplectic since we do not fully solve ([5]). In Fig.[5]below, we will observe 
that the no-loop version has comparable error behavior with a reduced computational cost, though it 
seems difficult to justify why the energy and invariants are well preserved in this version. 



12 



2.3 Numerical results 



In this section, we provide numerical tests of the behavior of our algorithm. As mentioned above, we 
are not interested in computing exact trajectories for the fast variables. Rather, we are interested 
in how well the algorithm preserves invariants of the system as well as its prediction of quantities 
derived from the fast variables. This system undergoes a slow exchange between the actions (|33l) of 
the fast degrees of freedom, and we test how well the algorithm captures this exchange. We also test 
the robustness of the algorithm by testing for the appearance of resonances. 

We compare our algorithm to a well-known integrator for highly oscillatory dynamical systems, 
the Mollify algorithm [7], which is a modification of the previously proposed Impulse (also known as 
Verlet-I/r-RESPA) algorithm [15] ■ Both Impulse and Mollify follow a kick/oscillate/kick pattern 
and incorporate the slow forces VU only at the "kick" steps, which are separated by a macro time step 
that is large with respect to the shortest timescale in the solution. This time step is typically larger 
than e and is thus larger than the stable regime for Verlet. For the "oscillate" step, these methods 
integrate the fast forces using a stepsize that is small with respect to e. The Mollify algorithm differs 
from the Impulse algorithm in how it incorporates the forces at the "kick" steps in order to improve 
the stability of the algorithm in the face of resonances, an issue we discuss later. These algorithms 
are designed to minimize the number of evaluations of the slow force, with the assumption that the 
"oscillate" step is cheaper, or comparable in cost, to a single evaluation of the slow force. In our tests, 
we have used the Verlet scheme for the "oscillate" step within Mollify, with inner stepsize equal to 
e/100. 

2.3.1 Modified FPU 

The Fermi-Pasta-Ulam (FPU) chain is a commonly-used test case for highly oscillatory integrators [TUl 
Sec. XIII. 2.1]. The chain is a collection of alternating stiff, harmonic springs and soft, nonlinear springs 
(see Fig. (JTJ) ) . After a change of coordinates, the potential can be written as in ([6]), with £1 constant. 




Figure 1: Fermi-Pasta-Ulam spring chain. 

We choose / = s = 3, corresponding to 3 stiff springs and 4 soft springs (there is one less degree 
of freedom since the total chain length is prescribed). For positions q\ = (91,1,91,2,91,3) € R 3 and 
92 = (92,1,92,2,92,3) € R 3 , the potential V is given by 



V(qi,q 2 ) = ^ f (91,1 ~ 92,i) 4 + ^(91,1+1 - 92,-i+i - 9m - 92,i) 4 + (91,3 + fa.a) 4 J • 

We modify the FPU problem to have non- constant fast frequencies. We consider a Hamiltonian of the 
form ©, choose 



0(9! ) = ^1 + 91,1, 

and use V above for the slow potential. 

The behavior of the modified FPU is qualitatively similar to the original. The fast and slow 
variables have timescales of 0(e) and O(l), respectively. In addition there is a slow exchange of the 
fast actions 

1 



Ij 2fi(fc) 



2 0(gi) 2 2 

P2,3 ~2 ?2,j 



j = 1,2,3, (33) 
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over long time periods of 0(e 1 ). The quantity 

/ = h + h + h (34) 

is an adiabatic invariant (see [TJ[2] an d QUI Theorem XIII. 6. 3]): it is an 0(1) quantity that is almost 
preserved over long times, with oscillations only of magnitude 0(e). 
For all numerical experiments, the initial conditions are 

g! = (1,0, Of, p! = (l,0,0) T , q 2 = (e,0,0) T , p 2 = (l,0,0) T . 
For these values, we have H £ (0) = 2.5 + 3e 2 + 0.5e 4 w 2.5 and 7(0) = - ( -= + V2 J w 1.06. 



2.3.2 Preservation of invariants and exchange of actions 

We first monitor how the energy (j9|) and the adiabatic invariant ()34|) are preserved along the numerical 
trajectory. Fig. [2] shows the lack of drift in both quantities over long times, for h — 0.005, e = 10~ 3 
and the above choice of initial conditions. Note that energy preservation for symplectic schemes is 
typically proven in the limit h — > 0, for fixed e. In our case, h/e — 5, so that energy preservation is 
not guaranteed by the theory. 



Energy and Adiabatic Invariant 

1.11 
1.1 

1.09 
1.08 
1.07 
1.06 
1.05 
1.04 



Be -1.4 

Adiabatic Invariant 



250000 500000 750000 le+006 
Long-time Trajectory 



Figure 2: Energy (we plot H £ — 1.4 to allow better scaling) and adiabatic invariant computed with 
Algorithm [TJ over the long time interval [0, 10 6 ] for e — 10~ 3 . We use stepsize h = 0.005. 



In Fig. [3J we examine the slow exchange among the actions Ij defined in (1331) . We observe that 
Algorithm [T] accurately reproduces exchange of the actions, in contrast to Mollify. Note that we 
have used a smaller stepsize h for Mollify than for Algorithm [T] in order to balance one possible 
measure of computational cost - the number of evaluations of the slow force. We discuss the issue of 
computational efficiency in Section 12.3.31 

2.3.3 Resonance and computational efficiency 

One common difficulty encountered among integrators for highly oscillatory systems is the appearance 
of resonances which destroy the long-time preservation of energy and other invariants. In linear cases, 
these typically occur when the slow time step h is an integer multiple of half the period of the fast 
motion. In Fig. 01 we explore the resonance behavior by simulating many trajectories, with various 
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Aig.ni 



Verlet 



Mollify 




2500 5000 7500 10000 
time 




2500 5000 7500 10000 
time 




2500 5000 7500 10000 
time 



Figure 3: A single trajectory is simulated using three algorithms, Algorithm [TJ Verlet, and Mollify. 
The energy, adiabatic invariant I = Ii + I2 + I3, and individual actions Ij are plotted to time 10 4 , 
with parameter e = 10~ 3 . Verlet is run with a very small time step h — 10~ 5 and is considered as 
the reference solution. Mollify and Algorithm [T] are run with stepsizes h — 0.0008 and h = 0.02, 
respectively, which are chosen so that both algorithms involve roughly the same number of calls to 
the slow forces. This plot shows that Algorithm Q] better captures the slow exchange of actions for a 
comparable computational effort. 



choices of the ratio h/e. We do this in two different ways: by holding e fixed while varying h as well 
as by holding h fixed while varying e. 

We display the maximum error in energy and maximum variance in / to time T = 10 4 , 

crr= max \H e (t) - H e (0)\, var = max lift) -1(0)1. (35) 
te[o,io 4 ] te[040 4 ] 

Recall that H e (0) w 2.5 and 1(0) « 1.06. Since the exact trajectory preserves energy, we desire that 
the method has small error in H 6 . On the other hand, there is variation in I(t) even for the exact 
trajectory and correctly predicting the variation is also of interest. Therefore, for each e, we have also 
computed a reference variation in I. Results are shown in Fig. 01 

For fixed s = 10~ 3 (the upper row of Fig.0]), Algorithm U performs quite well in the whole range 
of h, with only a few spikes in the error. Mollify only performs well for small h, with large, generalized 
resonant regions. For fixed h = 0.02 (the lower row of Fig. [4]), Algorithm [T] performs better for smaller 
e, and the error in the energy blows up for increasing e. The prediction of the variation in I is even 
better, matching the reference result. Note that the fact that Algorithm [T] performs better for smaller 
e is consistent with the fact that it is derived using homogenization techniques which are valid in the 
limit e — > 0. For large e, the terms neglected in the generating function are no longer small. 

We next consider the maximum error in energy versus the computational cost. For the sake of 
our comparison here, we use the number of evaluations of the slow force — VV as a proxy for overall 
computational cost. We thus assume that the slow forces are much more expensive to compute than 
the fast forces. For example, in a simulation of molecular chains, if the fast forces represent the bonds 
between adjacent particles in the chain and the slow forces include all other intramolecular as well 
as long-range intermolecular forces, then the slow force is much more computationally involved. For 
the Mollify algorithm, this will mean assuming that the computational expense of the 'oscillate' step 
is negligible in comparison to evaluating the slow forces, which may, depending on the application, 
underestimate the computational cost. The 'oscillate' step involves propagating the position and 
momentum according to the fast forces using a small time step and, in addition, propagating the 

Qq Qp Q £ p ar tial derivatives of the position and momentum with respect to the initial 



matrix 



P P 
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Error in Energy for e = 10 



Variation in I for e = 10 





h/e 



Error in Energy for h = 0.02 



Variation in / for h = 0.02 





Figure 4: Comparison of resonances for Algorithm [T] and Mollify. In the upper row, a series of 
trajectories are simulated for different stepsizes h, with e — 10 -3 . In the lower row, different values 
of e are considered, and the dynamics is integrated using the stepsize h = 0.02. In both cases, the 
trajectory is simulated to time T = 10 4 . In the left (resp. right) column, the maximum error in the 
energy (resp. the maximum variation in I) over the trajectory, as defined in (|35|) . is plotted. For the 
variation of /, a reference value is calculated using the Verlet algorithm with a very small time step. 
Here, Mollify exhibits many more resonances than it does in the case of constant f2, as shown in |1Q[ 
Chap. XIII] and |T2| Figs. 4-10]. 
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Alg. [TJ no loop 
Mollify 




-06 lc+07 le+08 le+09 
Gradient Calls 




-06 lc+07 le+08 
Gradient Calls 



le+09 



Figure 5: Efficiency plots comparing Algorithm [TJ a no-loop variant, and Mollify. For each graph, a 
series of trajectories are simulated for various stepsizes h with a single set of initial conditions and a 
single choice of parameter e (e — 10 -3 for the left plot and e = 10 -4 for the right plot). The trajectory 
is simulated to time T = 10 , and the maximum variation in the energy over the trajectory is plotted. 



condition [TJ [T3]. This matrix of size (s + /) x (s + /) is propagated according the Hessian of the 
fast potential in the so-called variational equation. In practice, simulating the variational equation 
is much more expensive than simulating a Hamiltonian dynamics on q 6 M. s+ f according to the fast 
forces. In the sequel, the cost of the 'oscillate' step in Mollify is entirely neglected. 

Figure [5] displays the maximum error in energy to time T — 10 4 versus the computational cost. 
Two plots are shown: we consider both cases e — 10 -3 and e = 10~ 4 , and integrate the dynamics with 
different stepsizes h. As the stepsize decreases, the computational cost increases. For all algorithms, 
the error is generally decreasing as a function of computational cost, up to the presence of resonances. 
In addition to Algorithm [TJ and Mollify, we also consider a variant of Algorithm [TJ that does not fully 
loop in order to solve the implicit equations (15Tj) (see (|3"2"jl ). 

Although Algorithm [TJ is much more stable with respect to resonances (as shown in Fig. [TJ , the 
trade-off is that the computational cost per time-step is increased. Thus, for e = 10~ 3 , Mollify is 
slightly cheaper than Algorithm [TJ in terms of cost required to achieve a certain accuracy. However, 
the computational cost of the no-loop variant is much smaller than that of Algorithm [TJ and is also 
smaller than that of Mollify. In addition, the no-loop variant does not exhibit significantly more 
resonances than does Algorithm [TJ For e = 10~ 4 , Algorithm [TJ enjoys a slight edge in comparison to 
Mollify, and the no- loop variant offers a significant advantage in computational cost. 

We also observe from the comparison of the cases e = 10~ 3 and e = 10~ 4 that, when e decreases 
and h is kept fixed, the computational cost of Mollify remains the same, while its accuracy decreases. 
On the other hand, the cost of Algorithm [TJ decreases somewhat since the number of iterations to 
solve the implicit equations (pTi"j) decreases, and the accuracy increases because the terms neglected in 
the expansion become smaller. In addition, resonances seem to disappear. 

3 The case of a matrix-valued, constant frequency 

In this section, we consider the case where the fast frequency is a constant matrix. We consider the 
Hamiltonian (see ([7])) 

rj , \ Pi Pi , P2P2 ., rl , . q%n 2 q 2 

H s [qi,q2,Pi,P2) = — — + —rj- + V{qi,q2) H , (36) 
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and assume that Q is a diagonal, positive definite constant matrix of size / X / (if f2 is a constant 
positive definite symmetric matrix, then, upon diagonalizing J7, we recover the case considered here). 
Our aim is to show that, upon a slight modification, we can also apply our strategy, thereby extending 
our previous study in |12) which was restricted to the scalar case Q = wld. In the sequel, we have to 
distinguish two cases, whether or not the eigenvalues of f2 are resonant. 

In both cases, we follow the exact same strategy as that we used in [12] and in Section [2) we 
first precondition the fast motion using a change of variables and next apply our two-scale ansatz to 
the Hamilton- Jacobi equation associated to the new Hamiltonian. The difference between the non- 
resonant case and the resonant case lies in the ansatz we make. In the non-resonant case, we introduce 
a fast time for each frequency in the generating function (see (|42|) V In contrast, in the resonant case, 
we introduce a unique fast time for each group of frequencies that are resonant one with each other 
(see |5B]) ). In both cases, the identification process follows the same lines. The derivation of the scheme 
is presented in Sections 13.11 and 13.31 in the non-resonant case and in the resonant case, respectively. 
Numerical results illustrating the non-resonant case are reported in Section 13.21 We conclude this 
section by considering a test-case with three distinct frequencies, two of which are resonant with each 
other (see Section |3~4|) . 



3.1 The non-resonant case: derivation of the scheme 

Without loss of generality, we assume that the matrix in (|36|) reads 



I a>ild mi 
w 2 Id TO2 



\ 



V 











(37) 



Wdld md / 



with < oj\ < . . . < LOd- The multiplicity of uji is m„ with y^-, m, = /. In view of the block diagonal 



decomposition of ft, we decompose qi £ W into d vectors 32,1 G 
and 

d 

Tr\2 \ 2 T 

q 2 il q 2 = 2^uJi g 2 ,i92,i- 



such that q 2 = (92,1, • ■ • , 92,, 



Likewise, we write p 2 = (p2,i> ■ • • ,P2,d)- We assume in the following that the tOi are non-resonant, in 
the sense that, for any k E Z d , 



LOjkj = => k = 0. 

3=1 



(38) 



As pointed out above, we follow here the exact same strategy as that we used in [12], or in Section[2] 
we first precondition the fast motion, and next consider the Hamilton- Jacobi form of the equations. 
We proceed first with the time-dependent change of variables (q2,i,P2,i) G K 2mi H> {x 2 ^,y 2)i ) = 
Xi(t,<l2,i,P2,i) G R 2m ' defined by 



<12. 



e 



X2. 



£ . ( LJit 

— sin 

uii \ e 



2/2. 



P2,', 



■ sin 



which reads, in a more compact form, 



q 2 = cos 



x 2 + efi 1 sin 



2/2, P2 



0J{t 

e 
fit 

E 



X2A + COS 



OJjt 

e 



2/2/ 



x 2 + COS 



e J 



2/2, 



(39) 



with x 2 = (#2,1) ■ • ■ , x 2 ,d) and y 2 = (2/2,1, • ■ • , 2/2, d)- The dynamics on (x 2 , 2/2) reads 

'fit 



x 2 = eil 1 sin 



d 2 V 



gi(i),cos I — I x 2 (t) + sil sin I — I y 2 



(*) 



y 2 = - cos 



nt\ 



e J 



)d 2 V 



qi (t),cos 



X2 (t) + efl 1 sin 



2/2 (t) 
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and the dynamics on [q±, x 2 ,pi, y 2 ) is a Hamiltonian dynamics with the time-dependent Hamiltonian 



/■/■noil , , 1 Pi P 1 1 TT/non-rcso /" j^jj ., ,. ,1 , i m 

H e {t,qi,x 2 ,pi,y 2 ) h vv e [—,..., —,qi,x 2 ,y 2 ) , (-Hi) 



where 

wr n - reso (n,...,r (2 , ft , a;2! y 2 ) = y 



gi, (cosri)x 2 ,i H (sinri)j/ 2 ,i, . . . , (cosT d )x 2 .d H (smr d )2/ 2 d 



We take the Hamiltonian as a starting point for our manipulations. In the sequel, we proceed 
with the construction of an approximation S e (h, qi,x 2 , Pi,Y 2 ) of the solution S e (h, qi,x 2 , P%,Y 2 ) to 
the Hamilton- Jacobi equation associated to (14TJ1) . for small times h. Observing that the variable x 2 is 
of order e, we first perform a change of variables and of unknown function: 

= ^x 2 and S £ (t,q 1 ,r 2 ,Pi,Y 2 ) = S e (t, q 1 ,eCl~ 1 r 2 , Pi,Y 2 ) , 




^- rW ^-* cs °(^,..., L ^, qi +d Pl S e ,en- 1 r 2 + d Y2 S ei Y 2 j , 4| , 

S e (0,qi,r 2 ,P 1 ,Y 2 )=0. 
We make the ansatz 

S e (t,qi,r 2 ,Px,Y 2 ) = S (t, r x , . . . , r d , q 1} r 2 , P 1; Y 2 ) + e5 x (i, n, . . . , r d , q lt r 2 , P 1: Y 2 ) (42) 
+ higher order terms in e k , k > 2, 

where the fast times r, are defined by 

r t = ^, l<i<d, (43) 

£ 

and where the functions (Sk)k>o are assumed to be 2-7T periodic with respect to each Tj. We also 
assume that the functions are of class C 1+a with respect to each Tj, with a € N, a > d. 

Remark 3.1. In the case where V does not depend on q 2 , the solution to (|4"TT) can be analytically 
identified and is indeed of the form (|42|). 

We now insert (|4"2")l in (|4*Tj) , identify the first ci variables of If™"-™ with the fast times {T~i} 1<i<d , 

and expand in powers of e. Based on ([5|), we have that X 2 = x 2 + dy 2 S e . The fast position X 2 is 
of order e, so So does not depend on Y 2 . The equation of order e" 1 in the expansion of (|4*T1) then 
becomes 



y^^j s n 5 = o. 



Using Lemma 13711 below, we deduce that 

yi<i<d, d n S =0. 

Lemma 3.1. Let /(t 1; . . . ,r d ) be a function that is 2tt periodic with respect to each Tj, and of class 
(ji+a w j±}i respect to each of its argument, with a £ N, a > d. Assume that 



3=1 



for a constant c, and a d-tuple {<^j}i<cj<d sa ^ s fy^ n 9 the non-resonance condition (J35J). Then the 
function f is a constant and c = 0. 
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Proof. As / is periodic, we can write it as its Fourier series: denoting r = (n, . . . , T4), we have 

/(t)= /feexp(ifc-r). 
Oi,...,fc ci )ez< i 

The assumption (l44l) reads 

i ^ ujfkkj exp(«fc ■ r) = c. 

3=1 (kx,-,k d )^ d 

As / is of class C 1+Q , we have that \fk\ < C(l + |fc|) _1 ~ Q for some constant C, where we recall 
\k\ = Ylj=i Hence, we have that \fkkj\ < 00, and the above sum is well-defined. We deduce 

that 

d 

Vfc E Z d , fc ^ 0, f k Wjfcj = 0, (45) 
3=1 

whereas the identification for k = leads to c = 0. We infer from (|4"5j) and assumption ([551) that 
./fc = for any fc 6 Z d , fc 7^ 0. So /(r) = /o- This concludes the proof. □ 

The equation of order e° reads 

d t S + 5>i0 n Si = + V(<Zi + d Pl S , 0). (46) 

Since So does not depend on r and Si is 27r periodic with respect to each Tj, we can again apply 
Lemma \3. II We thus infer from (l46l) that 



ftS = + V(gi + d Pl S , 0) (47) 

and d Ti S\ — for all 1 < i < d. Equation (|47[) is supplied with the initial condition So(t = 
0, <?i , J"2 , Pi ) = 0. For each r2, we thus recognize the Hamilton- Jacobi equation for the Hamiltonian 
function 

H 1 (q 1 ,p 1 ) = ^ + V(q 1 ,0). (48) 
So So does not depend on r<i. In the sequel, we will approximate So(t, qi, Pi) by 

S SE (t,gi,Pi) = S o (0,«i,Pi) + tftSo(0,gi,i\) = * + ^(?i,0) > ) , (49) 



which amounts to integrating the Hamiltonian dynamics generated by (|48f with the symplectic Eulcr 
algorithm. We have S (t) = S^ E (t) + 0(t 2 ). 

The sequel of the identification is not difficult. The bottom line is again as follows: since the {t{\ 
are independent variables in each Sk, we can integrate with respect to each one to split the equations 
and close the hierarchy. Following arguments similar to those presented in [TJ], we find that 

S 1 =0 (50) 
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and that S 2 (t) = Sf E (i) + 0{t 2 ), with 

d 
d 

+ V — (V a ,iV(?i + Wi,0)) T ((sin Ti)r 2>i - (cosTi)F a ,i) 

2 — 1 2 

< d 1 

2 — 1 L 

where the derivatives of V are, unless otherwise mentioned, evaluated at (<?i,0), and where V| jV is 
the Hessian matrix of V with respect to q 2 j. 

Observe that, in (|51l) . there is no term coupling components associated to different frequencies uti 
and uij, j 7^ i. This is reminiscent of the fact that, in the ansatz, the fast times Ti are independent 
variables, and that each Sk is 2n periodic with respect to each Ti. 

Consider now the approximation S £ (h) ~ 5'" on_reso (/i), with 

S™ n - ieso (h) := S* E (h) +eS 1 (h) +e 2 Sl E (h), 

where Sq E , Si and 5| E are respectively defined by (|49|) . ((50)) and (|5T|) . Using this approxima- 
tion, we obtain a symplectic algorithm in the variables (qi, x 2 ,pi, y 2 ). Returning to the original 
variables (q±, q 2 ,Pi,P2), we obtain the symplectic Algorithm [2] outlined below, which we denote by 

(Qi,q 2 ,pi,p 2 ) = *r n_rosol (9i.92,Pi,P2). 
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Algorithm 2 (Preconditioned Symplectic Scheme TCSol (qi,q2,pi,P2))- Set (qi,q 2 ,pi,p 2 ) = 
(<Zi > Q2 jPuP?)) r i = LUih/e and perform the following steps: 


1. Change of variables: set X2 — q 2 , 1/2 — P2, T 2 — Hx 2 / £. 


2. Solve for (Pi,Y 2 ] 


in the equations 




2/2 


d d 

= ^ + ^£-V 2 iV(qi,0)r2,i +eir S -^V2, l V(q 1 + hP x ,Q), 




Pi 


a 1 

= Pi + fcVi%,0) + e 2 V — VwV(ax,G)Y 2 ,i 

1— 1 * 


< 




1 

+ £ 2 V— ViaiVfei +/iPi,0)((siiiTi)>"2,i - (coSTi)y 2)i ) 

i=i w * 
d 1 

~ te V— Vi2^(<7i,0)V 2 ,i^((7i,0) 

2—1 * 






i, £ 2 d 1 

+ — j- E — (^V 1 2i2i^(gi,0)r 2 ,i + 3^iViaiai^(gi,0)y 2 ,i). 


3. Set Qi = 


= 9i + 


ftPi + fe 2 V — Vi 2i V(gi + ftPi.O) ( (sin Ti)r 2 ,j - (cos Ti)y 2)i ) • 

z— 1 * 


4. Set 






X 2 = x 2 + 


d 

^ 2 E 

i=l 


— V a ,iV(gi, 0) + fe 2 V —jV^Cgi, 0)r a ,< - e 2 V ^V 2 , a V(<Zi + Wi, 0). 


5. Return to the original variables: set r = Q/i/e and 




Q 


2 = (cos t)X 2 + 1 (sin r)y 2) P 2 = -- (sin t)X 2 + (cost) Y 2 . 


Set (<tf +1 ,<£+ 




= (Qi,Q 2 ,Pi,P 2 ). 



Note that, at step 2, we need to solve a system implicit in Z = (Pi,Y 2 ), which reads z = Z + 
hF(Z) + eG(Z) with z = (pi, y 2 ). In practice, we use a fixed point method, which converges in only 
a few iterations (three iterations in the test-case considered below). 

Remark 3.2. High-order derivatives of V appear in Algorithm^ As shown in Section |2"T21 above, they 
can be replaced by a finite difference approximation in the generating function. 

Neglecting all terms of order e 3 , the scheme ^ on ^ reso (qi,q 2 ,pi,p 2 ) is first order in h. A simple, 
well-known, manner to get a scheme of higher order is to consider the symmetric form 

I r\ r\ f ) / > \ iTfiion— reso2 / „ _ \ (.t.iioii- resol\ ,T f noii- resol / _ _ \ 

{Qi,Q 2 ,Pi,P 2 ) = w h {qi,q 2 ,pi,p2) = [v h / 2 ) v h/2 {qi,q 2 ,pi,P2), 

where *f>* denotes the adjoint of ^. This scheme, denoted Algorithm [3] in the sequel, is symplectic, 
symmetric and, neglecting all terms of order e 3 , of order 2 in h. 



22 



Algorithm 3 (Preconditioned Symplectic Scheme *^ on rcso2 (<7i, q2,Pi,P2)) ■ Set (qi,q 2 ,pi,p 2 ) 
(li , tfziPiiPT) an d perform the following steps: 

1. Set (^.Sa.Pi.Pa) = *^- re ~ 1 (gi, fc.Pi.j*). 

2. Set (Q 1 ,Q 2 ,P 1 ,P 2 ) = (^r r6So1 )* (Qi,Q 2 ,Pi,P2)- 
Set («r +1 .?2 +1 ,Pi +1 ,i« +1 ) = (Qi,Q 2 ,Pi,P 2 ). 



3.2 The non-resonant case: numerical results 

We consider a Hamiltonian of the form (|3"o| . with gi 6 I and 52 = (92.1, 92,2, 92,3) £ K 3 , and where 
the slow potential energy is 

92) = (c+ 92,i + 92,2 + 7<72,3) 4 + ^qhli + \q\ 

with c = 1 and 7 = 2.5. We choose f2 = diag (l, 1, y2j as the matrix of fast frequencies. 
Let Ij denote the energy associated to each fast degree of freedom: 

2 2 2 

_ P2,j ^ 9 2 ,j 

^- 2 + 2£ 2 . 1<J<*> 

with wi = ZU2 = 1 and ZU3 = \/2- We first note that here we use Ij to denote the fast energy, as 
opposed to Section [2] where it denotes the fast action (energy divided by frequency). Of course, for 
the constant frequency case here, these two quantities differ only by a multiplicative constant. In 
addition, we use here a different convention than that of Section 13.11 on how the eigenvalues of fi are 
numbered. It is well-known (see [TUJ Sec. XIII. 9]) that the quantities 

3 

I = ^2lj and h (52) 

2=1 

are adiabatic invariants of the dynamics. In contrast, I\ and I 2 , associated to the same frequency, are 
not adiabatic invariants, although their sum, I\ + I 2 =1 — 1$, is. 

We first choose e = 1/70 and h = lOe, and we monitor the evolution of the energy and adiabatic 
invariants up to time T = 10 6 on the numerical trajectory computed with Algorithm [3] Results are 
shown in Fig. [5J We observe no drift. 

For the same parameters, we show in Fig. [7] the evolution of Ij over the time window [0, 50]. As 
expected, I3 is preserved, as well as I\ + I 2 . We observe that Algorithm [3] correctly reproduces the 
exchange between I\ and I 2 . 

We now focus on the robustness of Algorithm [3] as e decreases. We set the time step to h — 0.02 
and consider the variations of the energy and of the adiabatic invariants (|52[) , 

max m) ~* m , max M max MtzlM, (53) 
te[o,io*] H(0) ' telcio 4 ] 1(0) ' te[o,io*] Z 3 (0) V ' 

over the time interval t £ [0,10 4 ], for stiffness e varying between 10~ 3 to 1. Results are shown on 
Fig. [8] Only a few resonances can be seen, and they are extremely peaked. In addition, the algorithm 
is very stable as e decreases to 0. 
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250000 500000 750000 le+06 
t 

Figure 6: Energy and adiabatic invariants I and I3 (for convenience, we plot H E — 2, I and J3 + 0.3) along 
the trajectory computed with Algorithm [3] (e = 1/70 and h — lOe). 




Figure 7: Preservation of the adiabatic invariants I1 + I2 and ^3, and exchange between I\ and I2, for e = 1/70 
(Algorithm [3] has been used with h = lOe). 



0.0075 



0.005 - 



0.0025 - 




0.1 



5 10 15 20 

h/e 



0.02 h 




1 


1 1 \ 










\ 

\ 
\ 






, % 
\ ^ 




- / — -^H^: 




1 1 1 









10 

h/e 



15 20 



Figure 8: Maximum variations (153|) of the energy (left) and of the adiabatic invariants I and 73 (right) on 
the time interval [0, 10 4 ], for several e (h — 0.02), for Algorithm [3] 
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3.3 The resonant case: derivation of the scheme 

In this section, we again consider the Hamiltonian (|36l) . with a diagonal matrix f2 £ R^ X A In contrast 
to Section 13.11 we assume here that some entries of the diagonal of fi are resonant. To simplify the 
notation and the analysis, we assume that is of the form (1371) with just d = 2, namely 



n = 



UJbldm b 



with m a + TOf, = /, < u> a < ojf,, and 



Lib = -Wo, with a £ N*, (3 £N*, q^/J, and a and (3 coprime. (54) 
a 

Likewise, we write q 2 = (q 2 , a ,q2,b) with q 2%a £ K ma and q 2i b £ W nb , and p 2 = (j)2, a ,P2,b)- 

As in Section [3~Tl we first consider the time-dependent change of variables (f39|) . which we write 



{x 2 ,y 2 ) — xitjQiiP?) S K ' . The dynamics on (qi,x 2 ,pi,y 2 ) is a Hamiltonian dynamics with the 
time-dependent Hamiltonian 

Hl es °(t,q u x 2 , Pl ,y 2 ) = P l?± + wr o (^,q u x 2 ,y 2 ) , 

2 \ ae I 



where 

Wr°(T, qi ,x 27 y 2 ) = V 



£ £ 

qi, (cos ar)x 2 , a H (s'm aT)y 2 . a , (cos f3r)x 2t b H (sinpY)?/^;, 

W a CJf, 



We let S £ (t, qi, x 2 , Pi, Y 2 ) denote the solution to the Hamilton- Jacobi equation associated with HI 
and perform the change of variables and of unknown function 

r 2 = ^x 2 and S e (t, q u r 2 , Pi, Y 2 ) = S e (t, q 1: rfTVa, Pi, Y 2 ) . 



The function SV satisfies 



dt S £ = tlSL + w** so ( w °t 



qi + d Pl S s ,en- l r 2 + d Y2 S E ,Y 2 ^j , 



2 V«e 7 (55) 

S e (0,qi,r 2 ,P u Y 2 ) = 0. 

We make the ansatz 

5 e (t,gi,r 2 ,Pi,y 2 ) = ^ (i,r,gi,r 2 ,P 1 ,y 2 )+e5 1 (i,r, g i,r 2 ,Pi,r 2 ) (56) 
+ higher order terms in £ k , k > 2, 

where the fast time r is defined by 

tuj a tuj b 

T = = "p-i ( 57 ) 

ae pe 

and where the functions (Sk)k>o are assumed to be 27r periodic with respect to r. 



Remark 3.3. In the case where V does not depend on q 2 , the solution to (|55l) can be analytically 
identified and is indeed of the form (l56l). 



We now insert ([55)) in (1551) . identify the first variable of W^ cso with the fast time r, and expand in 
powers of e. As in Section [3~T1 So is independent from Y 2 and r. The equation of order e° reads 

a t s + — a T Si = + v( qi + d Pl s , o). (58) 

a 2 
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Since So does not depend on r and Si is 2tt periodic with respect to r, we infer from (|58[) that 

dtS = + 7(gi + d Pl So, 0) (59) 

and d T Si = 0. Equation (|59l) is supplied with the initial condition So(£ = 0, q%, r 2 , Pi) — 0. We again 
recognize the Hamilton- Jacobi equation for the Hamiltonian function (|48|). In the sequel, So is thus 
again approximated by (|49[) . 



Following the same arguments as in |12j . we next proceed with the sequel of the identification, 
using the fact that, since a ^ j3 in ([51)1 , 

cos ar cos/3r dr = 0. 

As a consequence, at least at the orders of the expansion that we consider, no coupling appears 
between the terms associated to the frequency u a and those associated with the frequency Wb- We 
find that 

Si = (60) 

and that S 2 (t) = Sf E (t) + 0(t 2 ), with 

Sl B (t,T,qi,r 2 ,Pi,Y 2 ) = \{V 2 . a V) T Y 2 , a + \{V 2 , b V) T Y 2 , b 



\{^2. a V{qi +tP 1 ,0)) T ((sin ar)r 2 . a - {cosar)Y 2<a ) 



— (V 2 . 6 y(qi +iPi,0)) T ((sin /3r)r 2)6 - (cos /3r)y 2 , b ) (61) 

^ f^(v 2 . Q ^) T v 2 . a ^ + ^(v 2 , b y) T v 2 . h ^ 

^2 KaVL^2,a + C V 2,a^2,a) 



K fc v^yr 2ib + r 2 T b v^t/y 2 , fc ) 

6 



4w 2 



where the derivatives of V are evaluated at (qi,0) unless otherwise mentioned, and where V| a V is 
the Hessian matrix of V with respect to g 2jQ . 



Observe that, in (|61l) . there is no term coupling components associated to different frequencies. 
This is reminiscent of the fact that, at the 0(e 2 ) term that we consider here, coupling terms can only 
come from terms containing products of the form cos ar cos /3t, the average of which vanishes. In 
contrast, to identify S3, we need to handle terms of the form (cosar) s (cos/3r) r (with r + s = 3), 
whose average may not vanish (e.g. in the case a = 1, j3 — 2, s = 2 and r = 1). 

We next observe that the generating functions Sq E , Si and S 2 E that we have identified in this 
resonant case, defined by (j49|) . (f60|) and ([6"Tj) . where the fast time r is given by (|57p . are equal to 
the generating functions Sq E , Si and S 2 E identified in the non-resonant case, see (|49|) . ([50)) and (|51[) . 
where the fast times Tj are given by ()43j) . As a consequence, even though the frequencies u a and 
are here resonant, we again obtain the algorithms [2] and [3] proposed above. 

3.4 The resonant case: numerical results 



In Sections 13.11 and 13.31 we have derived algorithms to integrate the dynamics in the non-resonant 
case and in the resonant case, respectively. As explained at the end of the previous section, when 
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expanding up to 0(e 2 ), these two cases turn out to yield the same algorithm. A natural intuition, 
which is confirmed by a careful identification similar to the ones performed above, is that this algorithm 
is also appropriate in the case when the matrix of fast frequencies is of the form {57} , where the 
again satisfy < uj\ < LJ2 < ■ ■ ■ < ^d, but where, in contrast to the situation considered in Sections 13. II 
and 13.31 they are now possibly resonant, and d is arbitrary. 

Consider a system of the form (f36|) . with qi g R and 92 = (92,1, 92,2, 92,3, 92,4) € K , and where the 
slow potential energy is 

V(q%, 92) = (c + 924 + 92,2 + 92,3 + 792,4) 4 + ^9i92,i + ^9? 

with c = 1 and 7 = 2.5. We choose the matrix of fast frequencies to be f2 = diag (l, 1, V2, 2). A very 
similar test-case has already been studied in [3] and [TUl Sec. XIII. 9.1], where the chosen value of 
c = 1/20 causes the exchange between I\ and I2 to occur on a slower time scale. It obviously enters 
the framework of the current section, but not the framework of Section 13.11 as the first and the last 
frequencies are resonant. In this case, the adiabatic invariants are 

4 

I = Y, I J and J 3> ( 62 ) 

i=i 

where Ij is the energy associated to the jth fast degree of freedom: 

pi i wf 9n i 

with ZUi = Z02 — 1, UJ3 — V2 and UJ4 = 2. 

We first choose e = 1/70 and h = lOe, and we monitor the evolution of the energy and adiabatic 
invariants up to time T = 10 6 on the numerical trajectory computed with Algorithm [3] Results are 
very similar to those shown on Fig. [6j we do not observe any drift. 

We now focus on the robustness of Algorithm [3J as e decreases. We set the time step to h = 0.02, 
and consider the variations ([55)1 of the energy and of the adiabatic invariants (p2"]l . over the time 
interval t e [0, 10 4 ], for stiffness e varying between 10 -3 to 1. Results are shown on Fig. [HI When e 
decreases to 0, the algorithm performs better and better, except for a few peaked resonances. 




5 10 15 20 5 10 15 20 

h/e h/e 



Figure 9: Maximum variations (|53|l of the energy (left) and of the adiabatic invariants / and ^3 (right) on the 
time interval [0, 10 4 ], for several e (h — 0.02), for Algorithm [3] There is an increase in resonances compared 
to Fig. [8] though as before they are very tightly peaked. 
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We conclude this section by discussing the exchange of fast energies. In the case at hand, it is well- 
known (see [TU1 Sec. XIII. 9]) that some exchange occurs between I\, I2 (which are both associated to 
the same frequency uJ±) and I4 (which is associated to a frequency ZD 4 resonant with TDi). The exchange 
between I\ and I2 is 0(1) over timescales of 0(e^ 1 ), which is stronger than the exchange between I4 
and I\ + 12, which is only 0(e) over timescales of 0(e _1 ). For some choices of parameters, there is an 
O(l) exchange between I4 and 1\ +I2 that occurs on the long time scale 0(e~ 2 ). It turns out that, on 
the numerical trajectory computed using Algorithm [3j the exchange between I\ + 12 and I4 on long 
time scales is not reproduced, and I4 is almost preserved. This is due to the fact that (i) resonant 
frequencies are handled by the algorithm we derived as if they were non-resonant (see Section I3.3[) 
and (ii) in the non-resonant case, no exchange occurs between the fast energies. If we were to expand 
the generating function to third order in e, certain exchange terms among the fast terms would be 
included in the resulting algorithm, which could potentially capture the exchange between I\ + 12 and 
/;• 

Yet, the fact that we chose not to include these terms did not destroy the preservation of energy and 
adiabatic invariants, which are accurately recovered, as pointed out above. In addition, the exchange 
between I\ and I2, that occurs on the time scale 0(e _1 ), is accurately reproduced by Algorithm [3j as 
shown on Fig. QUI 
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Figure 10: The exchange between 7i and I2 on the time scale 0(e _1 ) is accurately reproduced by Algorithm[3] 
as well as the preservation of the adiabatic invariant I\ + I2 + Ii- The 0(e) exchange between I4 and Ii + I2 
is not captured. (We take e = 1/70, and we use h = lOe for Algorithm [3]) . 



4 The extensible pendulum 

In the previous sections, we have considered Hamiltonians of the form ©, where the leading order 
behavior of the fast variables is that of a harmonic oscillator, with either a slowly varying scalar 
frequency (in Section [2]), or a constant matrix of fast frequencies (in Section |3]). On the other hand, 
many examples of interest are not harmonic, and this technique cannot handle completely general fast 
forces. An intermediate class of system is the one discussed in [101 Sec. XIV. 3], where the Hamiltonian 
reads ^ ^ 

H(q,p) = -p T p + V a i ow (q) + U fast (q) 
with fast potential £/f as t that has a change of variables y — (2/1,2/2) = x(l) so that 

U iast (q) = ^V2^(yi) 2 V2, 

where ri(yi) is a / x / matrix. So, up to a change of variables, the fast potential energy is of the form 
considered in @. However, the map (q,p) h-> (y,p) is in general not symplectic. Hence we cannot 
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work with the variables (y,p), as the dynamics in these variables is not Hamiltonian. We thus need to 
also change variables for the momenta, so that the map (q,p) > (y, v), where v are the new momenta, 
is symplectic, and the dynamics in (y,v) is Hamiltonian. A standard consequence of this change of 
variables for the momenta is that the kinetic energy turns out to depend on the positions y. The 
Hamiltonian reads 

H(y,v) = T(y,v) + V slow ( X ^ (v)) + ^a^G/OV (63) 

Although the fast potential is harmonic, the fast Hamiltonian is not necessarily close to that of a 
harmonic oscillator. 

In this section, we consider a particular case (see (IM|) below) that cannot be written in the form ([5} . 
After a global change of variables, in positions and momenta, we transform the Hamiltonian to the 
form (|6"3")) (see below), where the kinetic energy turns out to still be a quadratic form of the 

momenta: T(y,v) = —v T M~ 1 (y)v. The resulting mass matrix M(y) depends on the position. We 
will show that our strategy can still handle this case and yields an efficient algorithm. 

4.1 Transforming to internal coordinates 

We consider an extensible pendulum in two dimensions (see Fig. The potential energy of the 

system is the sum of two terms, the spring energy and a term W(a) that depends on the angle a. 
Denote by q x and q y the Euclidean coordinates of the particle and by p x and p y the corresponding 
momenta. The Hamiltonian of the system reads 

H C3 , Ttesi3 , n {q x ,qy,p x ,Py) = ^{pl+p 2 y ) + -^(r - r ) 2 + W{a), (64) 

where r = ^Jq x + q y is the length of the spring, ro is the equilibrium length of the spring, e is a small 
parameter, and W is a 2-7r-periodic non-negative function. 

For illustration, this example can be considered as an oversimplified model of a molecular system, 
in which the spring-like potential (2e 2 ) _1 (r — ro) 2 models each bond length between neighboring 
atoms, and with some potential associated with each bond angle. 



(<lx,qy) 

Figure 11: Extensible pendulum test-case: a particle at position (q x , q y ) is attached to the origin using 
an extensible spring. We denote by r the spring length and by a the angle between the pendulum and 
the vertical. 

The fast oscillating term (2e 2 ) (r — ro) 2 of the Hamiltonian (|64|) is not harmonic with respect 
to the cartesian degrees of freedom (q x ,Qy), and thus the strategy described in Section [2] cannot 
be directly applied to (jo"4"]l . However, as is often the case in molecular dynamics models, the fast 
oscillating term of (|64j) is harmonic in the internal degrees of freedom, the bond length r and the 
bond angle a. It is thus natural to introduce polar coordinates (a,r), with 

q x =r cos a and q y —r sin a, 

so that the fast oscillating term of the potential energy is harmonic in these coordinates. We now 
introduce momenta (p a ,Pr) associated to (a,r), such that the map (q Xl q y ,p x ,p y ) H> (a,r,p a ,p r ) is 
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symplectic. Following jTOl Example VI. 5. 2], this leads to choosing 

Pr = Px cos a + p y sin a and p a = —p x r sin a + p y r cos a. 
In addition, rather than working with r, we work in the sequel with 

r = r -r , 

so that the fast oscillating harmonic term attains its minimum at r = 0, as in Without loss of 
generality, we furthermore assume that ro = 1. The Hamiltonian then reads 

H(a,r,p a ,p r ) = | + + ^r 2 + W(a). (65) 

Although the mass matrix depends on the fast position r, we will show that our general strategy, 
originally developed for can still be used. 

Remark 4.1. As above, we choose initial conditions depending on e such that the energy is bounded. 
As W is bounded from below, this implies that, at any time t, r(t) = 0(e). As a consequence, one 
could approximate (|65p by 



2 2 j 

fldeo D up(o,r,p a> p P ) = y + y + + W(a) 



In this case, (r,p r ) and (a,p a ) are decoupled. In the following, we do not make this approximation, 
and we work with (1651) . 



4.2 Derivation of the symplectic scheme 



Starting from (|65[) . we now follow our usual strategy: we first precondition the fast variables, and 
next consider the Hamilton- Jacobi form of the equations. Observe that, in the case at hand, the fast 
frequency is constant. We thus follow the approach described in [12] and consider the time-dependent 
change of variable (r,p r ) h-> (b,Pb) defined by 

y t . t 1 . t t 

r = ocos — hepf,sm- and p r = — sm — hp&cos-. (66) 

e e e e e 

In the variables (a,b,p a ,pb), the dynamics is again a Hamiltonian dynamics for the Hamiltonian 
# P cnd f~, a, b,p a ,Pb^j , with 

H pend (T,a,b,p a ,p b ) = — = : — + W(a). 

2(l + o cos r + epb sm t) z 

Let S e {t,a,b,P a ,Pb) solve the Hamilton- Jacobi equation associated to -Spend- As r and b are of 
order e, we make the change of variables and of unknown function: 

b=- and S £ (t,a,b, P a , Pb) = S E (t,a,eb,P a ,P b ) , 

so that S s satisfies 

d t S e =H peild (^,a + dp a S s ,eb + d Pb S £ ,P a ,P b \ 5 e (0, a, 6, P OJ P b ) = 0. (67) 



We make the ansatz 



S e (t,a,b,P a ,P b ) = S {t 1 T,a,b,P a ,Pb)+sS 1 {t,T,a,b,P ai P b ) (68) 
+ higher order terms in e k , k > 2, 



30 



where the fast time r is defined by 



t 



and where the functions (Sk)k>o are assumed to be 2n periodic in r. We now insert in (joTf . 
identify the first variable of -Spend with the fast time r, and expand in powers of e. 

The identification follows the same lines as in the previous sections. We obtain that So is inde- 
pendent of b, Pb and r, and satisfies So(t) — So(t) + 0(t 2 ) with 



S$ E (t,a,P a )=t 



P 2 
2 



W(a) 



while 



Si =0. 

Using the method of characteristics, we obtain that S2(t) = Sf E (t) + 0(t 2 ), with 
Sl E (t,T,a,b,P a ,P b ) = P 2 [P b cosT-bsmT}- P b (P a +tW'(a)) 2 
+ t 



\{b 2 + P 2 )(P a + tW'(a)) 2 l -(P a + tW'{a)f 



(69) 
(70) 
(71) 



Collecting (|69l) . (I7D1) and (I7T|) . we obtain the following approximation of the generating function: 



S e (h, a, b, P a , P b ) « S c (h, a, b, P a , P b ) = S^(h, a, P a ) + e 2 S^(h, -,a,b, P a , P b 

Returning to the variables (a,r,p a ,p r ) of the Hamiltonian in (|65|) . we obtain a symplectic scheme, 
called Algorithm@]in the sequel, which we denote by (i 



,n+l r n+l „n+l ~,n+l 



,P n a + \p n r +1 ) = * 



pendulum/ n n 



Algorithm 4 (Symplectic scheme \]/P endulum ( a) r,p a ,p r )). Set (a,r,p a ,p r ) = (a™ 
and perform the following steps: 


,r n ,p%,p?), t = h/e 


1. Change of variables: set b — rje. p b =p r . 




2. Solve for (P a ,P b ) in the equations 






Pb = Pb- 


eP 2 sinr + \eb (P a + hW'{a)) 2 , 




< 


Pa = Pa + 


hW'(a) - 2he 2 P b {P a + hW'{a))W"{a) 






+ h 2 e 2 

- 


^{b 2 + P 2 )(P a + hW'(a)) - 2(P a + hW'(a)f 


W"(a). 


3. Set 










A = a + hP a + 2e 2 P a {P b cosT-bsmT)-2e 2 P b {P a + hW'(a)) 




+ he 2 


^(b 2 + P 2 )(P a + hW'(a)) - 2{P a + hW'{a)f 




4. Set B = b + 


eP 2 cosr- e{P a + hW'(a)) 2 + ~heP b (P a + hW'(a)) 2 . 




5. Return to the original variables: R — eB cos r + ePb sin r, P r = —B sin 


T + P b COS T. 


Set (a n+1 ,r n+1 


n n+l n+l\ _ 
i-Fq i-Fr J 


(A,R,P a ,P r ). 





To solve the implicit equations for (P a ,P b ), we first observe that, once P a is known, P b can be 
explicitly determined. We thus first recast the problem as a scalar implicit equation on P a and solve 
it using a fixed point algorithm before determining P b . 



31 



Neglecting all terms of order e 3 , the scheme ij/P ondulum j s g rs ^ order in h. As in Section ETTl we use 
pendulum to build a symplectic and symmetric scheme, denoted Algorithm [5] in the sequel, following 

t A r> r> r> \ lTr pend.sym. / \ /, T ,pcndulum\ * . T .pcndulum/ \ 

{A,R,P a ,P r ) = (a,r,p a ,p r ) = (#£ /2 J *f l/2 (a,r,p a ,p r ), 

where ^P endulum ^ denotes the the adjoint method, which is denoted Algorithm [Bl 



Algorithm 5 (Symmetric Symplectic Scheme \j/P ondsym ( aj r,p a ,p r )). Set (a,r,p a ,p r ) 
and perform the following steps: 



(a",r",K,K) 



1. Set (A,R,P a ,P r ) = *^ mnm (a,r,p a ,p r ). 

2. Set (A, R, P a ,P r ) = fe C /2 dUlUm )* (A,R,P a ,P r 



Set (a n+1 ,r n+1 - n+1 - n+1 



Algorithm 6 (Symplectic Scheme fijrP eIldulum j (a, r,p a ,p r )). Set (a,r,p a ,p r ) — (a™ , r n , p™ , p") , r = 
ft,/e and perform the following steps: 

r . r 

1. Rotate the variables: set b = - cos r + p r sinr, p& = — sin r + p r cos r. 

£ £ 

2. Solve for (A, P) in the equations 
' a = A - hp a + 2e 2 p a (p b cosT + Psinr) - 2e 2 p b (p a - hW'(A)) 



he 2 



\{B 2 +p 2 ){p a - hW'(A)) - 2( Pa - hW'(A)f 



2 

b = B + sp 2 a cos T-s{ Pa - hW'(A)) 2 - ^-hEp b {p a - hW'(A)) 2 . 



3. Set P b =p b + ep 2 a sinr - -heB{p a - hW'(A)) 2 . 



4. Set 



P a = Pa-hW'(A) + 2he 2 p b ( Pa ~hW'{A))W"(A) 
"3, 



2_2 



h 2 e 



^■{B 2 +p 2 )( Pa - hW'(A)) - 2(p a - hW\A)f 



W"{A). 



5. Return to the original variables: R — eB, P r = P b . 
Set (a n+1 ) r n+1 ,P2 +1 ,P? +1 ) = (A,R,P ai P r ). 



4.3 Numerical results 

We have implemented Algorithm [3] on the Hamiltonian (|65p with 

W» = (cos a) 2 . 

We first choose e — 2.10 -3 and h = 0.02, and monitor the evolution of the energy and of the adiabatic 
invariant 



1 9 

2 2e 2 



(72) 
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up to time T = 10 6 . Results are shown on Fig. Q21 We observe no drift. 



We next study the robustness of the algorithm as e decreases. We set the time step to h = 0.02, 
and consider the variations (|53[) of the energy and of the adiabatic invariant (|72[) . over the time 
interval t E [0, 10 4 ], for stiffness e varying between 10 -3 to 1. Results are shown in Fig. [121 Again, 
the algorithm performs equally well when e decreases to 0, up to extremely peaked resonances in the 
energy preservation. 
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Figure 12: Left: energy (for convenience, we plot H — 0.9) and adiabatic invariant I along the trajectory 
computed with Algorithm [5] (e = 2.10 and h = 0.02). Center and right: maximum variations (|53|l of the 
energy (center) and of the adiabatic invariant / (right) on the time interval [0, 10 4 ], for several e (h = 0.02), 
for Algorithm [5] 
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